// Copyright 2020 Poofee (https://github.com/Poofee)
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// ------------------------------------------------------------------------
/*****************************************************************************
 *                                                                           *
 *  File:    relay1250.cpp                                                   *
 *                                                                           *
 *****************************************************************************
 *                                                                           *
 *  Authors: Poofee                                                          *
 *  Email:   poofee@qq.com                                                   *
 *  Address: https://github.com/Poofee                                       *
 *  Date:    2020年07月15日                                                   *
 *                                                                           *
 *****************************************************************************/
#include "relay1250.h"

#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#if !defined(ARMA_32BIT_WORD)
#define ARMA_32BIT_WORD
#endif
#include <armadillo>
#include <vector>
#include <ctime>
#include <omp.h>

#include "SuperLU_MT.h"
#include "slu_mt_ddefs.h"


using namespace std;
using namespace arma;


Relay1250::Relay1250() {

}

Relay1250::~Relay1250() {

}


/*!
 * @brief      参数初始化
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay1250::init() {
	/** 文件名，需不带扩展名 **/
	setFileName((char*)"../geo/JRS1250/JRS1250bgmband");	

	MAX_NONLINEARSTEPS = 200;
	/** 设置材料参数 **/
	CMaterial* air_material = new CMaterial;
	materialList.emplace_back(air_material);

	CMaterial* pm_material = new CMaterial;
	materialList.emplace_back(pm_material);
    pm_material->H_c = 0;//1.3 / 1.12 / (4 * PI * 1e-7);
    pm_material->miu *= 5000;//1.12;

	CMaterial* dt4e_material = new CMaterial;
	materialList.emplace_back(dt4e_material);
    dt4e_material->sigma = 0;//1 / (2e-7);

	int num_bh = 63;
	dt4e_material->BHpoints = num_bh;
	dt4e_material->miu = 3000 * miu0;
	dt4e_material->Hdata = (double*)malloc(num_bh * sizeof(double));
	dt4e_material->Bdata = (double*)malloc(num_bh * sizeof(double));
	/** 1250的BH曲线 **/
	double H[] = { 0,10,100,150,200,250,300,350,360,370,380,390,400,410,420,430,440,450,460,470,480,490,500,510,520,530,540,550,560,570,580,590,600,610,620,630,640,650,660,670,680,690,700,710,720,730,740,750,760,770,780,790,800,900,1900,2900,3900,4900,5900,6900,7900,8900,9900 };
	double B[] = { 0,0.0599718475592522,0.597090304698867,0.891647197655098,1.17932987869717,1.44849641137579,1.66011089327454,1.76457710473790,1.77545758128628,1.78436476327606,1.79172378239890,1.79786552945965,1.80304417672406,1.80745433989896,1.81124544595357,1.81453294220234,1.81740674973144,1.81993754230050,1.82218138737627,1.82418318028023,1.82597919696382,1.82759900417761,1.82906689986371,1.83040300831993,1.83162411995517,1.83274434065704,1.83377559810748,1.83472803974086,1.83561034796116,1.83642999167750,1.83719342844891,1.83790626803611,1.83857340558121,1.83919913072030,1.83978721749903,1.84034099887956,1.84086342880472,1.84135713415656,1.84182445846260,1.84226749882734,1.84268813727457,1.84308806745612,1.84346881750189,1.84383176964217,1.84417817711895,1.84450917881082,1.84482581192233,1.84512902302871,1.84541967771808,1.84569856903387,1.84596642488713,1.84622391458191,1.84647165457468,1.84851935249073,1.85578788972131,1.85843034336099,1.86030242999031,1.86190677974238,1.86338688482529,1.86479924261073,1.86617060520324,1.86751527801887,1.86884160416511 };
	for (int ih = 0; ih < num_bh; ++ih) {
		dt4e_material->Hdata[ih] = H[ih];
		dt4e_material->Bdata[ih] = B[ih];
	}

	CMaterial* coil_material = new CMaterial;
    materialList.emplace_back(coil_material);
    Coil* coil = new Coil;
    coil->xyz       = 1;
    coil->direction = 1;
    coil->Nc        = 60;
    coil->Rc        = 3.9;
    coil->Vc        = 27;
    coil->type      = CICULARCOIL;
    coil->height    = 17e-3;
    coil->center[0] = 0;
    coil->center[1] = 0;
    coil->center[2] = 0;
    coil->r0        = 0;
    coil->width     = 3e-3;
    /** 线圈磁场应向下-27 **/
    coil->Jor = coil->Vc / coil->Rc * coil->Nc / coil->width / coil->height;
    coil_material->coil = coil;
    coil_material->tau = coil->Nc / coil->width / coil->height;

	/** 设置映射关系 **/
	setMaterial(1, dt4e_material);/** 衔铁 **/
	setMaterial(2, coil_material);/** 线圈 **/
	setMaterial(3, dt4e_material);/** 铁芯 **/
	setMaterial(4, dt4e_material);/** 铁芯 **/
	setMaterial(5, dt4e_material);/** 铁芯 **/
	setMaterial(6, pm_material);/** 永磁 **/
	setMaterial(7, dt4e_material);/** 铁芯（导磁环） **/
	setMaterial(8, dt4e_material);/** 铁芯 **/
	setMaterial(9, air_material);/** 外部空气 **/
	setMaterial(10, air_material);/** 空气 **/
	setMaterial(11, air_material);/** 空气 **/
	setMaterial(12, air_material);/** 空气 **/
	setMaterial(13, air_material);/** 可压缩空气 **/
	/** 设置形变区域 **/
	setXiantieTag(1);
	setAirTag(9);/** 空气开始的编号，要保证后面都是空气 **/
	setRemeshTag(13);/** 最后一个区域重分网 **/

	Precision = 1e-6;
    ParallelThreads = 4;
}


/*!
 * @brief      使用传统的牛顿法进行求解瞬态模型
 * @author     Poofee
 * @date       2020年9月18日
 * @param[out] 
 * @return     void  
**/
void Relay1250::MagnetoDynamics2DAxisNR() {
	/** 计时 **/
    double t0, t1, t2;
	clock_t time[10];
	int tt = 0;
	time[++tt] = SuperLU_timer_();
	/** 电路参数 **/
    double C1 = 200e-6;
    double C2 = 200e-6;
	double R1 = 1e3;
	double Rc = 3.9;
    double Y1 = 1 / R1;
    double Yc = 1 / Rc;

	double ce[3][3] = { {0} };/** 雅可比矩阵 **/
	double Je[3][3] = { {0} };/** 用于牛顿迭代 **/
	double Te[3][3] = { {0} };/** 涡流矩阵 **/
	double De[3]    = {  0  };/** 电流密度矩阵 **/

	double Y11, Y12, Y13, Y22, Y23, Y33;
	/** 变量 **/
    //double beta = 1;
	/** 时间步循环 **/
	current_step = 0;
	double error;
	double totalTime = 0;
    int timesize = 2;//timesteps.size();
    for (int i = 0; i < timesize; ++i) {
        printCurrentTime();
		double dt = timesteps[i];
		totalTime += dt;
		totaltimes[i] = totalTime;
        printf("Time step %d/%d, current time is %lf s.\n", i, timesize, totalTime);

        /** remesh，要使用增量位移，向下为负，向上为正 **/
		remesh(i);

		/** 读取第i步的分网 **/
		loadMesh();
        unitConvert(unitRatio);

		/** 查找边界，默认一次边界 **/
		findBoundaryEdges2D(-1);
        findBoundaryPoints2D(-1);
        if(0) {
            testBoundaryResults((char*)"../geo/JRS1250/bnd.vtk");
        }
		
		A_last = A;
        A_old = A;

		allocateMemory();
		setBoundary();
		findAir();
		initMu();

		/** 自由度：未知节点+电流 **/
		int num_dof = allPoints.size() - boundaryPoints.size();
        printf("DoF %d\n", num_dof+1);
		double ucoil = 0;//
        if(i > 0) ucoil = UCoil[i - 1];
		double phicoil = 0;
		double* unknown_b = (double*)calloc(num_dof+1, sizeof(double));
		double miuold = 1;
		double c = ((Y1 + Yc) * dt + (C1 + C2)) / (2 * PI);

		for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
			int i1 = num_ele - num_triangle + i_tri;
			/** 计算除了磁导率部分的系数 **/
			makeTriangle(i1);
			/** 计算感应电流的一部分D^{T} **/
			CElement* ele = pmeshele + i1;
			CMaterial* mat = getMaterial(ele->geometry_tag);
            if(mat->coil) {
                double dtmp = mat->tau * ele->AREA / 3.0;
                for (int j = 0; j < 3; ++j) {
                    bbJz(ele->n[j]) += dtmp;
                }
            }

		}
		bool LinearFlag = true;
		/** 非线性迭代求解 **/
        SimulationTime[current_step] = 0;
		for (int non_iter = 0; non_iter < MAX_NONLINEARSTEPS; ++non_iter) {
            t2 = SuperLU_timer_();
            printCurrentTime();
            printf("Time step (%d/%d),\tNewton's step (%d/%d)\n", i, timesize, non_iter, MAX_NONLINEARSTEPS);
			/** 有限元装配 **/
			int pos = 0;
			bn.zeros();
			unknown_b[num_dof] = 0;
			phicoil = 0;
			for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
				int k, m, n;
				int i1 = num_ele - num_triangle + i_tri;
				CElement* ele = pmeshele + i1;
                k = ele->n[0]; m = ele->n[1]; n = ele->n[2];

				/** 计算除了磁导率部分的系数 **/
                makeTriangle(i1);
				
				double ydot = ele->ydot;
				double miut = ele->miut;/** 当前非线性迭代中的mu **/
				miuold = ele->miu;/** 上一时刻的mu **/
				/** 相当于C矩阵 **/
				Y11 = ele->Y11; Y22 = ele->Y22; Y33 = ele->Y33;
				Y12 = ele->Y12; Y13 = ele->Y13; Y23 = ele->Y23;

				/** 全局矩阵 **/
				ce[0][0] = Y11;      ce[1][1] = Y22;      ce[2][2] = Y33;
				ce[0][1] = Y12;      ce[0][2] = Y13;      ce[1][2] = Y23;
				ce[1][0] = ce[0][1]; ce[2][0] = ce[0][2]; ce[2][1] = ce[1][2];

				/** 计算雅可比矩阵 **/
				double v[3];
				v[0] = Y11 * A(k) + Y12 * A(m) + Y13 * A(n);
				v[1] = Y12 * A(k) + Y22 * A(m) + Y23 * A(n);
				v[2] = Y13 * A(k) + Y23 * A(m) + Y33 * A(n);

				CMaterial* mat = getMaterial(ele->geometry_tag);
				double tmp = 0;
				if (non_iter > 0 && mat->BHpoints > 0) {
					tmp = mat->getdvdB(ele->B);
					miut = mat->getMiu(ele->B);
					if (ele->B > 1e-9) {
						tmp /= ele->B * ele->AREA;  // B==0?
						tmp /= ydot * ydot * ydot;
					}
				}
				if (mat->BHpoints > 0) {
					LinearFlag = false;
				}
				Je[0][0] = v[0] * v[0] * tmp; Je[1][1] = v[1] * v[1] * tmp; Je[2][2] = v[2] * v[2] * tmp;
				Je[0][1] = v[0] * v[1] * tmp; Je[0][2] = v[0] * v[2] * tmp; Je[1][2] = v[1] * v[2] * tmp;
				Je[1][0] = Je[0][1];          Je[2][0] = Je[0][2];          Je[2][1] = Je[1][2];

				/** 计算右侧向量 **/
				double dtmp = 0;
				De[0] = dtmp; De[1] = dtmp; De[2] = dtmp;
                double dtmp1 =  -2 * PI * Yc * dtmp;

                if (mat->coil) {
                    dtmp = mat->coil->direction*mat->tau * ele->AREA / 3.0;
                    De[0] = dtmp; De[1] = dtmp; De[2] = dtmp;
                    dtmp1 = 2 * PI * Yc * dtmp;
					/** 计算线圈磁通 **/
                    phicoil += 2 * PI * dtmp * A(k);
                    phicoil += 2 * PI * dtmp * A(m);
                    phicoil += 2 * PI * dtmp * A(n);
                }

				/** 计算涡流矩阵 **/
                double etmp = mat->sigma / ydot * ele->AREA / 12;
				Te[0][0] = etmp * 2 + dtmp1 * bbJz(k); Te[1][1] = etmp * 2 + dtmp1 * bbJz(m); Te[2][2] = etmp * 2 + dtmp1 * bbJz(n);
				Te[0][1] = etmp * 1 + dtmp1 * bbJz(m); Te[0][2] = etmp * 1 + dtmp1 * bbJz(n); Te[1][2] = etmp * 1 + dtmp1 * bbJz(n);
				Te[1][0] = etmp * 1 + dtmp1 * bbJz(k); Te[2][0] = etmp * 1 + dtmp1 * bbJz(k); Te[2][1] = etmp * 1 + dtmp1 * bbJz(m);

				for (int kk = 0; kk < 3; ++kk) {
					for (int mm = 0; mm < 3; ++mm) {
						Je[kk][mm] += ce[kk][mm]/(miut * ydot);
                        if (i > 0) {
							Te[kk][mm] /= dt;
							Je[kk][mm] += Te[kk][mm];
                        }
					}

                    /** 上一非线性迭代步的值 **/
                    bn(ele->n[kk]) -= ce[kk][0] / (miut * ydot) * A(k) + ce[kk][1] / (miut * ydot) * A(m) + ce[kk][2] / (miut * ydot) * A(n);

                    if (i > 0) {
                        bn(ele->n[kk]) -= Te[kk][0] * A(k) + Te[kk][1] * A(m) + Te[kk][2] * A(n);
                        bn(ele->n[kk]) -= -Yc * De[kk] * ucoil;
                    }

                    /** 上一时刻的值 **/
                    if (i_tri < air_position && i > 0) {
                        bn(ele->n[kk]) += Te[kk][0] * A_last(k) + Te[kk][1] * A_last(m) + Te[kk][2] * A_last(n);
                    }

                    /** 计算永磁部分 **/
                    bn(ele->n[kk]) -= mat->H_c / 2. * ele->Q[kk];
                }

				/** 装配 **/
				int n_row, n_col;
				for (int row = 0; row < 3; ++row) {
					n_row = node_pos(ele->n[row]);
					for (int col = 0; col < 3; ++col) {
						/** 判断节点是否在未知节点内 **/
						/** 得到排序之后的编号 **/
						n_col = node_pos(ele->n[col]);
						if (n_row < num_dof && n_col < num_dof) {
                            insertSparseMatrixData(pos++, n_row, n_col, Je[row][col]);
						}
					}
                    if(i > 0){
                        insertSparseMatrixData(pos++, n_row, num_dof, -Yc * De[row]);
                        insertSparseMatrixData(pos++, num_dof, n_row, -Yc * De[row]);
                    }
				}
                unknown_b[num_dof] -= -Yc * De[0] * A(k) - Yc * De[1] * A(m) - Yc * De[2] * A(n);
				
                /** 主要是为了防止A_last访问越界 **/
                if (i_tri < air_position && i > 0) {
                    unknown_b[num_dof] += -Yc * De[0] * A_last(k) - Yc * De[1] * A_last(m) - Yc * De[2] * A_last(n);
				}
			}/** 单元循环结束 **/

			/** 右下角的点 **/
			if (i > 0) {
                insertSparseMatrixData(pos++, num_dof, num_dof, c);
			}
			if (non_iter == 0) {
				locs.reshape(2, pos);
				vals.reshape(1, pos);
			}

			/** 调用线性求解器求解 **/
			int size = num_dof;
			if (i > 0)
				size += 1;

            t0 = SuperLU_timer_();
			/** 使用构造函数来生成稀疏矩阵 **/
			sp_mat X(true, locs, vals, size, size, true, true);

			for (int idof = 0; idof < num_dof; ++idof) {
				unknown_b[idof] = bn(node_reorder(idof));
			}

			/** 对应电压的部分 **/
			if (i > 0)
                unknown_b[num_dof] += -c * ucoil + (C1 + C2) / 2 / PI * UCoil[i - 1];

            t1 = SuperLU_timer_();
            printf("Assembly time:%lf s\n", t1 - t0);

            if(0) {
                char fn[256];
                sprintf(fn, "%s%dunknown_b%d.txt", fileName, i, non_iter);
                FILE *fp1 = fopen(fn, "w");
                for (int it = 0; it <= num_dof; ++it) {
                    fprintf(fp1, "%lf\n", unknown_b[it]);
                }
                fclose(fp1);
            }
            printf("unknown_b[num_dof]:%10.8e\t\n", unknown_b[num_dof]);

            //---------------------superLU_MT---------------------------------------
            CSuperLU_MT superlumt(size, X, unknown_b, ParallelThreads);
			if (superlumt.solve() == 1) {
				printf("Error: superlumt.slove. Info:%d\n", superlumt.info);
				return;
			} else {
				A_old = A;
				double* sol = superlumt.getResult();

				for (int idof = 0; idof < num_dof; ++idof) {
					pmeshnode[node_reorder(idof)].A = sol[idof];
					A(node_reorder(idof)) = sol[idof];
				}
                error = norm((A), 2);
                if(error > 1e-10)
                    error = norm((A), 2) / norm(A_old, 2);
				A += A_old;
				if(i > 0)
                    ucoil += sol[num_dof];

                if(0) {
                    char fn[256];
                    sprintf(fn, "%dsol%d.txt", i, non_iter);
                    FILE *fp1 = fopen(fn, "w");
                    for (int it = 0; it <= num_dof; ++it) {
                        fprintf(fp1, "%lf\n", sol[it]);
                    }
                    fclose(fp1);
                }
			}

			/** 更新磁场结果 **/
			double bx, by;
			for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
				bx = 0; by = 0;
				int i1 = num_ele - num_triangle + i_tri;
				CElement* ele = pmeshele + i1;
				for (int j = 0; j < 3; ++j) {
					bx += ele->Q[j] * A(ele->n[j]);
					by -= ele->P[j] * A(ele->n[j]);
				}
				CMaterial* mat = getMaterial(ele->geometry_tag);
				ele->B = sqrt(bx * bx + by * by) / 2. / ele->AREA / ele->ydot;
				ele->Bx = bx / 2. / ele->AREA / ele->ydot;
				ele->By = by / 2. / ele->AREA / ele->ydot;
				ele->miut = mat->getMiu(ele->B);
                ele->miu = ele->miut;
			}

			/** 更新电路结果 **/
			if (i > 0) {
				//ICoil[i] = -(C1 + C2)*(ucoil - UCoil[i - 1]) / dt - ucoil * Y1;
				ICoil[i] = ucoil * Yc - Yc*(phicoil - PhiCoil[i - 1]) / dt;
				printf("E:%10.8e\n", (phicoil - PhiCoil[i - 1]) / dt);
			}				
            t0 = SuperLU_timer_();
            printf("Solving time:%lf s\n", t0 - t1);
			printf("error:%10.8e\tucoil:%10.8e\ticoil:%10.8e\n", error,ucoil, ICoil[i]);
			printf("Phicoi:%lf\n", phicoil);
            SimulationTime[current_step] += t0-t2;
			if (error < Precision || LinearFlag) {
                IterationStep[current_step] = non_iter;
				break;
			}
			bn.zeros();
			pos = 0;
		}/** 牛顿迭代结束 **/

		/** 计算电磁力 **/
        double forceY = calcMagForce(-1, FORCEMETHOD::MAXWELLTENSOR);
        /** 加速度->速度->位移 **/

		/** 计算线圈磁通 **/
		phicoil = 0;
		for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
			int k, m, n;
			int i1 = num_ele - num_triangle + i_tri;
			CElement* ele = pmeshele + i1;
			k = ele->n[0]; m = ele->n[1]; n = ele->n[2];

			/** 计算除了磁导率部分的系数 **/
			makeTriangle(i1);
			CMaterial* mat = getMaterial(ele->geometry_tag);
			if (mat->tau > 0) {
				phicoil += 2 * PI * mat->tau * A(k) * ele->AREA / 3;
				phicoil += 2 * PI * mat->tau * A(m) * ele->AREA / 3;
				phicoil += 2 * PI * mat->tau * A(n) * ele->AREA / 3;
			}
		}
		printf("Phicoi:%lf\n", phicoil);

		for (int i_tri = 0; i_tri < num_triangle; ++i_tri) {
			int i1 = num_ele - num_triangle + i_tri;
			pmeshele[i1].miu = pmeshele[i1].miut;
        }

		/** 导出数据 **/
        char fn[256];
        sprintf(fn, "%s%d%sAnoair.vtk", fileName, i, prefixName);
		writeVtkFileNoAir(fn);

        /** 更新结果 **/
		PhiCoil[current_step] = phicoil;
        magforces[current_step] = forceY;
		//Ddisplacements[current_step] = 0;
		//displacements[current_step]  = 0;
		if (i > 0) {
			velocities[current_step]     = Ddisplacements[current_step] / dt;
			accelerations[current_step]  = (velocities[current_step] - velocities[current_step-1]) / dt;
			
			UCoil[current_step]          = ucoil;
			//ICoil[current_step]          = (C1+C2)*(UCoil[current_step-1] - UCoil[current_step])/dt- UCoil[current_step]*Y1;
		}	
		printf("ICOIL:%lf\n", ICoil[current_step]);

		stepIncrement();
		
		/** 回收空间 **/
		if (unknown_b) delete unknown_b;
	}
}


/*!
* @brief      BBTLM求解二维瞬态
* @author     Poofee
* @date       2020年8月30日
* @param[out] 
* @return     void  
**/
void Relay1250::MagnetoDynamics2DAxisBBTLM() {

}


/*!
 * @brief      1250静态参数初始化程序。
 * @author     Poofee
 * @date       2020年10月25日
 * @param[out] 
 * @return     void  
**/
void Relay1250::StaticInitial() {
	printf("1250 relay static initial\n");
    init();
    setPrefixName((char*)"static");
}

/*!
 * @brief      1250瞬态参数初始化程序。
 * @author     Poofee
 * @date       2020年10月25日
 * @param[out] 
 * @return     void  
**/
void Relay1250::TransientInitial() {
	printf("1250 relay transient initial\n");
	init();
    setPrefixName((char*)"eddy");

	/** 初始化时间步长 **/
	current_step = 0;
	min_time = 0;
	max_time = 5e-3;
	timesteps.emplace_back(0);
	for (double t = 1e-5; t < 1e-3; t += 2e-5) {
		timesteps.emplace_back(2e-5);
	}
	for (double t = 1e-3 + 5e-5; t <= 5e-3; t += 5e-5) {
		timesteps.emplace_back(5e-5);
	}

	/** 初始化变量 **/
	current_step = 0;
	initPhysicalVariables(timesteps.size());
    UCoil[current_step] = 420;/** 第一步求解永磁 **/
	UCoil[current_step + 1] = 420;
	/** 初始化位移，用来测试和验证 **/
	double dis[] = {0,      0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2.57776e-09,-1.40519e-08,-4.01657e-08,-8.65132e-08,-1.58572e-07,-2.61644e-07,-4.00719e-07,-5.81379e-07,-8.09056e-07,-1.08895e-06,-1.42639e-06,-1.82780e-06,-2.30018e-06,-2.85057e-06,-3.48663e-06,-4.21614e-06,-5.04765e-06,-5.99038e-06,-7.05415e-06,-8.24941e-06,-9.58731e-06,-1.10797e-05,-1.27387e-05,-1.45772e-05,-1.66087e-05,-1.88474e-05,-2.13081e-05,-2.40062e-05,-2.69583e-05,-3.01810e-05,-3.36917e-05,-3.75083e-05,-4.90702e-05,-6.29622e-05,-7.95024e-05,-9.90193e-05,-0.000121834,-0.000148268,-0.000178617,-0.000213188,-0.000252266,-0.000296112,-0.000344961000000000,-0.000399048000000000,-0.000458567000000000,-0.000523709000000000,-0.000594631000000000,-0.000671459000000000,-0.000754324000000000,-0.000843317000000000,-0.000938524000000000,-0.001030,-0.001130,-0.001220,-0.001320,-0.001410,-0.001510,-0.001600,-0.001690,-0.001780,-0.001870,-0.001960,-0.002050,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100,-0.002100 };
	for (int i = 1; i < timesteps.size(); ++i) {
		displacements[i] = dis[i];
		Ddisplacements[i] = displacements[i] - displacements[i - 1];
	}
}

/*!
 * @brief      1250静态计算
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay1250::StaticAnalysis() {
	StaticInitial();
}

/*!
 * @brief      1250瞬态计算
 * @author     Poofee
 * @date       2020年10月24日
 * @param[out] 
 * @return     void  
**/
void Relay1250::TransientAnalysis() {
	TransientInitial();

	openGeo();

    MagnetoDynamics2DAxisNR();

	/** 输出信息 **/
	outputResults();
}

